Regularization of fields for self-force problems in curved spacetime: 
foundations and a time-domain application 
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We propose an approach for the calculation of self-forces, energy fluxes and waveforms arising 
from moving point charges in curved spacetimes. As opposed to mode-sum schemes that regularize 
the self-force derived from the singular retarded field, this approach regularizes the retarded field 
itself. The singular part of the retarded field is first analytically identified and removed, yielding a 
finite, differentiable remainder from which the self-force is easily calculated. This regular remainder 
solves a wave equation which enjoys the benefit of having a non-singular source. Solving this wave 
equation for the remainder completely avoids the calculation of the singular retarded field along 
with the attendant difficulties associated with numerically modeling a delta function source. From 
this differentiable remainder one may compute the self-force, the energy flux, and also a waveform 
which reflects the effects of the self-force. 

As a test of principle, we implement this method using a 4th-order (1-1-1) code, and calculate the 
self-force for the simple case of a scalar charge moving in a circular orbit around a Schwarzschild 
black hole. We achieve agreement with frequency-domain results to ~ 0.1% or better. 

PACS numbers: 04.25.D-, 04.25. dg, 04.25. Nx, 04.20. Cv. 04.30.Db 
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I. INTRODUCTION 

With the advent of gravitational wave astronomy ap- 
proaching, the development of accurate and efficient 
models for gravitational wave sources has steadily pro- 
gressed. The ability to predict gravitational wave ampli- 
tudes and waveforms for expected sources will greatly en- 
hance the usefulness of detectors such as VIRGO, LIGO, 
and LISA. An interesting class of relevant sources in- 
cludes a large 1O^-1O^°M0 black hole in a binary with a 
closely orbiting stellar-mass compact object. 

The orbit and inspiral of a compact object into a sub- 
stantially more massive black hole presents a complica- 
tion for traditional numerical analysis. A numerical grid 
must be fine enough to resolve the geometry in the vicin- 
ity of the small object, where the metric appears to be 
that of the compact object with tidal distortions from 
the large hole. But the grid must be coarse enough to 
reach the wave-zone of the binary, so that the waveforms 
might be carefully monitored. In addition, the timescale 
for the effect of radiation reaction is long compared with 
the orbital period and the waveform from many orbits 
will be used in the data analysis. Together the dramat- 
ically different length scales coupled with the dramati- 
cally different time scales present a formidable challenge 
for the study of an extreme mass ratio inspiral (EMRI). 

Perturbative analysis appears more feasible for the 
EMRI problem. The orbiting compact object is modeled 
as a point mass whose motion, to lowest-order, approxi- 
mately follows a geodesic in the background geometry of 
the large black hole companion. Thus, the emitted grav- 
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itational waves can be calculated reasonably well using 
the mature perturbation theory of black holes [U, 0, 0] . 
But in this approach it is necessary to go beyond the ap- 
proximation of geodesic motion in the background geom- 
etry. The effects of radiation reaction on the orbital phase 
requires an extension of the usual perturbation analysis 
to include what is often called the self-force. 

In a so-called frequency-domain approach, one chooses 
to Fourier decompose the source and the field and then 
solves for each Fourier mode of the field independently. 
This method works well for a flux calculation if the spec- 
trum is simple, such as that of a particle in a circular 
orbit. However for generic trajectories, including partic- 
ularly those which reflect the effects of radiation reaction, 
the frequency spectrum is complicated enough to make 
the frequency-domain analysis numerically expensive. 

Further, the field of the particle is singular at the par- 
ticle's location, and some regularization procedure is re- 
quired to calculate self-force effects. The usual procedure 
to date is termed mode-sum regularization, as initially de- 
veloped by Barack and Ori [3, H, 0] . This regularization 
prescription depends crucially upon a decomposition of 
the derivatives of the field into angular modes, such as 
spherical harmonics, which are individually finite. From 
each finite mode, the part that contributes to the sin- 
gularity but not the self-force is identified and removed; 
the remainders of the field derivatives for each mode are 
then summed to determine the finite effect of the par- 
ticle's field on its own motion. Generally the sum has 
convergence which is only polynomial in the mode num- 
ber and, thus requires analysis at high mode numbers 
0, i, i, IS inl, [13, m, H m to obtain accurate results. 

In this manuscript we introduce a general method for 
analyzing the field of a point charge orbiting a black hole 
and for directly determining the waveforms and flux in- 
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tegrals as well as the instantaneous self-force acting back 
on the charge itself, which includes all of the effects of 
radiation reaction in a natural manner. 

The strength of this approach lies in the derivation of a 
wave equation for a regular field -0^ which is identical to 
the retarded field in the wave zone and whose derivatives 
at the charge determine the self-force. We call the deter- 
mination oiip^, field regularization. Solving our effective 
wave equation requires neither Fourier decomposition in 
time, nor any angular decomposition for treating the dra- 
matically different length scales in the EMRI problem, 
and circumvents the need for ever calculating the actual 
singular retarded field. 

The source 5efF of the effective wave equation follows 
from a local analysis of the singular part ^p^ of the re- 
tarded field, Eq. Importantly, this effective source is 
smooth everywhere except for its limited differentiability 
at the location of the charge. Ample freedom in choosing 
S'off allows the source to spread out over a region with a 
length scale comparable to the size of the black hole or 
even to the distance from the charge to the black hole. 

We shall describe our approach in terms of a point 
source with a scalar charge interacting with its own 
scalar field while orbiting a large black hole. The for- 
mal extension of these ideas to pure gravity with a small 
Schwarzschild black hole perturbing the geometry of a 
much larger black hole is completely straightforward at 
the perturbative level. The details of this extension to 
gravity are algebraically complicated but conceptually 
simple and will be the focus of a future report. 



II. ORGANIZATION OF THIS PAPER 

The main objective of this work is to provide a proof of 
principle for the process of field regularization described 
in mill as a time-domain technique for self-force calcula- 
tion. It verifies that we are able to achieve results com- 
parable to that obtained with frequency-domain methods 
[9j, or other time-domain methods relying on the mode- 
sum decomposition [ll, • 

In ijlVI we describe the details of our numerical im- 
plementation of field regularization applied to a scalar 
charge in a circular orbit in Schwarzschild. Tests of the 
internal consistency of the numerical implementation are 
in ^ 

Section IVII displays the results of our self- force analy- 
sis for a scalar charge in circular orbits at Schwarzschild 
radius R = IQM and 12M. The time and radial compo- 
nents of the self force are compared with results from a 
frequency-domain analysis. We also reconstruct the en- 
tire retarded field and compare this with the retarded 
field of the frequency domain analysis. 

The discussion in Will summarizes our results and de- 
scribes the strengths of field regularization in comparison 
with other methods of self-force calculation and also with 
methods of current interest for calculating energy fluxes 
and waveforms for generic orbits about black holes. 



Appendix [XI gives some details of the expansion of the 
singular field ^jJ^ about the point charge and describes 
how higher order terms in the expansion increase the 
overall efficiency of field regularization. 



III. FIELD REGULARIZATION 

For a scalar charge, the general strategy for computing 
the self-force first involves solving the minimally-coupled 
scalar wave equation with a point charge q source. 



Airq / S'^'^\x - z{T))dT, 



(1) 



for the retarded field 'ip''°^. Here Va is the derivative 
operator associated with the metric gat of the background 
spacetime and 7 is the worldline of the charge defined 
by 2:°(r) and parameterized by the proper time r. The 
physical solution of the resulting wave equation will be a 
retarded field that is singular at the location of the point 
charge. As such, a self-force naively expressed as 



(2) 



will need a regularization prescription to make sense. 
Early regularization prescriptions [ll, [13, [11] were based 
upon a Hadamard expansion of the Green function, and 
showed that for a particle moving along a geodesic the 
self force could be described in terms of the particle in- 
teracting only with the "tail" part of ip, which is finite 
at the particle itself. Later [I^l it was realized that the 
singular part of the field which exerts no force on the 
particle itself could be identified as an actual solution to 
Eq. (IT|) in a neighborhood of the particle. A formal de- 
scription of in terms of parts of the retarded Green's 
function [T^ is possible, but generally there is no exact 
functional description for -0"^ in a neighborhood of the 
particle. Fortunately, an intuitively satisfying descrip- 
tion for tp^ results from a careful expansion about the 
location of the particle: 



^ps ^ qjp + Oipyn^) as p ^ 0, 



(3) 



where 7?, is a constant length scale of the background 
geometry and p is a scalar field which simply satisfies 
p2 —2:^ + 1/^-1- in a very special Minkowskii-like lo- 
cally inertial coordinate system centered on the particle, 



first described by Thorne, Hartle and Zhang 
applied to self- force problems in Refs. [1, [2 
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surprisingly the singular part of the field, which exerts 
no force on the particle itself, appears as approximately 
the Coulomb potential to a local observer moving with 
the particle. 

Our proposal for solving Eq. ([T]), and determining the 
self-force acting back on the particle now appears ele- 
mentary. First we define 



0^ = q/p 



(4) 
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as a specific approximation to tp^ . By construction, we 
know that tp^ is singular at the particle and is else- 
where. Also, within a neighborhood of the worldline of 

the particle 



as p — > 0. 



(5) 



Next, we introduce a window function W which is a 
C°° scalar field with 



= 1 + 0{p^l'R}) as p ^ 0, 



(6) 



and — > sufficiently far from the particle, in particular 
in the wavezone. Finally we define a regular remainder 
field 



which is a solution of 



(7) 



V"Va(VFVi^)-47r(7 / j('*)(x-z(T))dr (8) 



from Eq. H]). 

The effective source of this equation 



V''Va(VFVi^)-47rg / b''^\x - z{T))dT (9) 



•Scff : 



is straightforward to evaluate analytically, and the two 
terms on the right hand side have delta-function pieces 
that precisely cancel at the location of the charge, leaving 
a source which behaves as 



5cff = Oip/n^) as p^Q. 



(10) 



Thus the effective source 5cff is continuous but not neces- 
sarily differentiable, C", at the particle while being C°° 
elsewhere [S^]- Fig. [T] shows the source function which is 
actually used in the numerical analysis described in i jlVI 
The modest non-diffcrentiability of SeS at the particle is 
revealed in Fig. [21 
A solution of 



(11) 



is necessarily at the particle, and its derivative there 
provides the self force acting on the particle. Also, in the 
wavezone W effectively vanishes and ip^' is then identi- 
cally V''^"* and provides both the waveform as well as any 
desired flux measured at a large distance. 

General covariance dictates that the behavior of S'cff in 
Eq. ^ may be analyzed in any coordinate system. But, 
only in the specific coordinates of Refs. [20] and [2l| is it 
so easily shown Q that the simple expression for ip^ in 
Eq. dS]) leads to the O^pjV}) behavior in Eq. ^ and 
then to the nature of the solution of Eq. pip . 

We describe the procedure of solving Eq. (fTT|) as field 
regularization. Then the derivatives of "0^ determine the 



self-force, and is identical to ■0'"''' in the wave zone. 
With this process there is no apparent reason to deter- 
mine the actual retarded field. However, if one wants 
to compare results from field regularization with results 
from a traditional determination of the retarded field 
then simply adding Wip^ to the remainder results 
in the retarded field -0'^°*. Such a comparison for our 
trial of field regularization appears in Figs. (O and pU)). 




FIG. 1: The effective source Soft on the equatorial plane. The 
particle is at r/M = 10, (j)/n = 0, where Sett appears to have 
no structure on this scale. The smooth "double bump" shape 
far from the charge is a characteristic of any function similar 
to V^(VK/|r — ro|) in flat space, with a window function W as 
given in Eq. dTS]) . 




FIG. 2: The effective source Soft in the equatorial plane in the 
vicinity of the point source at r/M = 10, (f>/n = 0. Note the 
significant difference of scales with Fig. [1] 
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IV. NUMERICAL IMPLEMENTATION 



where V{r^) is implicitly given in terms of r as: 



As a concrete example and test of the field regu- 
larization prescription, we apply it to the well-studied 
case of a scalar charge moving in a circular orbit about 
a Schwarzschild black hole. We choose q/m — 1 for 
the charge to mass ratio of the particle, and a cir- 
cular geodesic at Schwarzschild radii, R = lOM and 
R = 12A'I, where M is the mass of the black hole. We 
work in Schwarzschild coordinates, in which the met- 
ric is expressed as gab — diag(— (1 — 2M/r), 1/(1 — 
2M/r), r2, sin^ 6*)'. 

For this task, we have developed code that (a) solves 
the regularized wave equation (jlip and (b) computes the 
scalar self-force. For simplicity, we have chosen to solve 
the regularized wave equation using a (l-l-l)-approach. 
We exploit the spherical symmetry of the background, 
decompose physical quantities into spherical harmonics, 
and then solve the resulting set of (l-l-l)D-wave equations 
(one 'time' -I- one 'space') for the spherical-harmonic 
components. 

It must be stressed at this point that the numerical 
implementation presented in this paper does not high- 
light the advantages of our prescription. The simplicity 
of the orbit wc consider and the spherical symmetry of 
our background geometry naturally lend themselves to 
a significantly more efficient frequency-domain approach. 
But the point here is to provide a quick, first check of our 
ideas. One is cautioned not to let the simplicity of the 
present problem obscure the generality of our proposed 
method, and its potential for cases with generic orbits 
and spacetimes lacking symmetry, and for self-consistent 
evolutions which are likely to require self-force calcula- 
tions in real time (as opposed to being a post-processing 
step). Testing the robustness of our method against these 
more difficult problems will be addressed in future work. 
The current goal is mainly to establish plausibility: to 
provide both an initial proof-of-principle for the method 
and also the necessary practice en route to tackling more 
interesting problems handled using more sophisticated 
numerical techniques. 



A. Scalar fields in a Scliwarzschild geometry 

Wave equations in spherically-symmetric backgrounds 
simplify considerably with a spherical-harmonic decom- 
position of the field. In the case of a Schwarzschild ge- 
ometry expressed in Schwarzschild coordinates, this de- 
composition is typically performed as follows: 



(12) 



With = r + 2A/ln(r/2M— 1), this yields equations 
for fimit, n): 



drl 



/^_2MX 



while the source Sim{t,r^,) is 



1(1 + 1) 



2M 



(14) 



Sim{t,r,) = {r-2M) / p{x'^)Yi^{6' ,,p')dn' . (15) 



In a frequency-domain approach, one further chooses to 
Fourier-decompose fimit,r^) = / F/mtj(r*) exp(— iwt)(ia;, 
and thereby solve the resulting set of ordinary differential 
equations for i^;mii;('"*)j for each mode lo. This method 
tends to be numerically expensive, however, for sources 
with a continuous a;-spectrum. Instead, we choose to 
solve Eq. (jl3p as an initial boundary value problem, in a 
time-domain fashion, for each {l,m). This is done with 
Sim computed beforehand as the spherical harmonic com- 
ponents of the effective source found in Eq. ([9]). 



B. Effective source term 

A novel feature of our approach is the use of an ef- 
fective source that permits the easy calculation of both 
self-forces and fluxes. As discussed above, this effective 
source is formally 

S^,, = ^V^{W-4!^) - Airq I 5^^\x - z{T))dT. (16) 
To lowest order, the singular field takes on the form 



(17) 



We take advantage of the results in where p is 
expressed explicitly as p = ^J'qtjX^x^ in Thorne-Hartle- 
Zhang coordinates for a particle moving in a circular or- 
bit. Using the coordinate transformation found in Ap- 
pendix B of [§], where a more detailed discussion of the 
singular field is found, we are able to express the singular 
field in Schwarzschild coordinates. (A brief discussion of 
this coordinate transformation is provided in Appendix 
A). To complete our effective source, we select a window 
function whose role is to kill off smoothly the singular 
field in regions where it is not needed. Consequently, the 
effective support of the windowed singular field Wij^^ is 
confined to a compact region surrounding the particle's 
world line. 

Our chosen window function is spherically-symmetric 
with respect to the center of the black hole. This choice 
was not necessary but guarantees that W would not un- 
necessarily modify the (/, m)-spectrum of the source, and 
thereby allows us to make more controlled comparisons 
with existing frequency-domain results on the same prob- 
lem. Our simple choice of W is 



V{r^)Jlra = Slra{t,r^) (13) 



W{r) 



exp 



R) 



N 



tN 



(18) 
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In this window function, the constant a sets the width, 
and the exponent N controls how quickly W and VqPF 
reach the required values of 1 and 0, respectively, as one 
approaches the particle. We use cr = 2M and N — 8 
in all the results presented in this paper. It is necessary 
that N is an even integer, and taking full advantage of 
the accuracy of our approximation for requires that 
W = 1 + 0(p^/7^'*) as p -> 0. Thus we require that N > 
4. In fact we used = 8 in anticipation of improving 
the approximation for in the future. 

Our choice for the window function leads to the effec- 
tive source ScS displayed in Figs. ([T]) and 0. A larger 
choice for a would spread the bumps out further, and a 
smaller choice for N would smooth the bumps. But if N 
were less than 4, then W^p^' would not adequately match 
the behavior of ip^ as p ^ 0. 

With the effective source constructed as above, its 
spherical-harmonic components were then computed. 
Circular orbits proved advantageous here because of 
which the time dependence of the components could then 
simply be inferred. The spherical harmonic components 
were evaluated with a 4th-order Runge-Kutta integrator 
with self-adjusting step size, which was derived from a 
routine in [24]. 



C. Evolution algorithm 

The integration scheme we use in evolving Eq. (jl3[) 
follows a technique first introduced by Lousto and Price 
[25| . and later improved to fourth-order accuracy by 
Lousto [1^ and Haas [13] ■ Unlike their schemes, how- 
ever, we do not deal with sourced and vacuum regions 
of our numerical domain separately. Their use of a sin- 
gular delta- function source meant that the resulting field 
was non-differentiable at the location of the charge, while 
smooth everywhere else. For us, the effective source is 
C*', implying that the field is at least C^. While this 
is still of finite differentiability, we find that the effective 
source is diffcrcntiable enough not to warrant a treatment 
different from the vacuum case. 

In the (t, r,)-plane, we introduce a staggered grid with 
step sizes At = ^Ar* — h. In this grid, a unit cell 
is defined to be the diamond region with corners {{t + 
h, r*), {t — h, r*), {t, r* -I- h), {t, r* — h)}. Only at these 
grid points do we evaluate We henceforth drop the 
spherical-harmonic indices in fim for convenience. 

The main idea behind the algorithm is to integrate the 
wave equation over a unit cell. This is done easiest with 
Eddington-Finkelstein null coordinates 7i = t — and 
V = t + as the integration variables. 

The differential operator of the wave equation, when 
expressed in {u,v) coordinates, is just — 4(9,j9„. Over a 
unit cell then, the derivative term in Eq. (|13p can be 





{t - ft.,r«) 

FIG. 3: Staggered (characteristic) grid with unit cell, 
integrated exactly: 




-4 dud, fdudv = -4[/(i + h, r,) + f{t - h, r,) 



- f{t,r,+h)- f{t,r,~h)]. (19) 

Integrations of the potential term and the source term 
do not enjoy the same simplicity as the derivative term. 
We need to approximate these integrals to the appro- 
priate order in h so as to achieve the desired 0{h'^)- 
convergence over the entire numerical domain. 

Suppose we wish to solve the wave equation over a 
region defined by AT and Ai?*. In this region, there will 
be iV = ATAR^,/h^ cells. Achieving 0(ft.'*)-convergence 
for evolution means that we need to integrate the wave 
equation with an over- all error of at most 0(/i^) over the 
entire computational domain. For a unit cell, this means 
an approximation with an error 0{h^)/N ^ 0{h^). 

Such an approximation is achieved with the double 
Simpson rule. Consider a sufficiently diffcrcntiable func- 
tion G(i, r*) to be integrated over a unit cell. The double 
Simpson rule then reads: 

jj Gdudv= [ao„c„ + WGit,r,) 

+ 4{G{t + h/2, - h/2) + G{t + h/2, + h/2) 
+ Git - h/2, - h/2) + G{t - h/2, + h/2))] 
+ 0{h'^), (20) 

where Gcomcn, is just the sum of the values of G evaluated 
at the corners of the unit cell. 

This is directly applied in integrating the source term 
of Eq. 

jj^Sf^dudv. (21) 
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One simply evaluates the source term at the required 
points and then sums these accordingly in order to get 
an 0(/i^)-accurate approximation to the integral. 
However, for integrating the potential term: 



-Vf du dv, 



(22) 



we recall that one has only restricted access to /. The 
direct evaluation of / is done only at the grid points, 
i.e. corners of the unit cell. Thus far, only G^^^^^^ in 
Eq. (PD|) can be explicitly evaluated. To use Eq. (^0)1 for 
the potential term, we need to determine how to evaluate 
/ at all the other points. 

Following Lousto [IHl, we evaluate G — —Vf at the 
central grid point (i.e. G{t, r*)) using values at the neigh- 
boring grid points on the same time slice. 

G(t,n) =^[9G{t,r, - h) + 9Git,r, + h) 

- G{t, n - 3h) - G{t, n + 3h)] + 0{h^). 

(23) 

Note that this is different from Haas [3], who uses grid 
points in the causal past of the unit cell. The 0(/i^ )-error 
incurred in this approximation is tolerable because of the 
/i^-factor that appears in Eq. (PU)) . 

We seek similar approximations for G in the remain- 
ing points. Consider first the pair G{t + h/2,r^, — h/2) 
and G{t — h/2,r^, — h/2). (The other pair, composed 
of G{t + h/2,r^ + h/2) and G{t ~ h/2,r^ + h/2), is 
treated similarly). This pair makes up the top and bot- 
tom corners of a smaller cell, Ci^ft , made up of the points 
{{t + h/2, - h/2), {t-h/2, - h/2), {t, - h), {t, r,)}. 

What we shall do next is find an approximation for 
G{t + h/2, - h/2) + G{t - h/2, - h/2) (24) 

accurate to 0{h^). Again, this is sufficient because of the 
/i^-factor in Eq. ([SO]) . 

Consider integrating the wave equation over this 
smaller cell, but this time only up to an accuracy of 
©(/i-*). The integral over the derivative term will again 
be exact: 



Cleft 



Adudyfdudv = -A[f{t + h/2,n - h/2) 
+ fit - h/2, - h/2) - fit, r,-h)- fit, 



(25) 

The integrals of the potential and source terms over this 
smaller cell are again handled as before, but this time 
we approximate them only to 0(/i*). To this end, the 
double trapezoidal rule will suffice, which reads: 



G du dv 



Cleft 



[Git + h/2,r^~h/2) 



+Git - h/2, - h/2) + Git, r^-h) + Git, r,)] 
+Oih^). (26) 



(*+ 1- •■•+!) 




-'right 



i) 



1) 



FIG. 4: Unit cell of the algorithm. The black dots indicate 
grid points, whereas the gray ones stand for the points where 
G — —Vf needs to be approximated. The subcells Cje/t and 
Cright are shaded gray. To approximate G at some of the gray 
dots, we integrate the wave equation in each of these subcells. 



Applying this to the potential term then gives: 
-Vfdudv = -('^) X 

Cleft V ^ / 

[y(r* - h/2) fit + h/2, - h/2) 
+ y(r* - h/2) fit - h/2, - h/2) 
+ Vir,-h)fit,r,-h) 
+ Vin)fit,r,)] + Oih^). (27) 

Combining Eq. ((25|) and Eq. ((27|) , the result of integrating 
the wave equation over this smaller cell yields: 



fit + h/2, n - h/2) + fit - h/2, n - h/2) = 

2 

ifit,r.,-h) + fit,r,)) 1 



Vir^ - h/2) 



5'eft dudv + Oih^). 



(28) 



Cleft 



After multiplying both sides of this last equation by 
— Vir„ — h/2), the resulting left-hand-side becomes two of 
the as yet missing pieces in the double Simpson formula: 

Git + h/2, - h/2) + Git - h/2, - h/2) = 
- Vir, - h/2) fit + h/2, - h/2) 
-Vir^-h/2)fit-h/2,r^-h/2). (29) 

The resulting equation then gives us the desired Oih^)- 
approximation of the missing expression. Git + h/2, — 
h/2) + Git - h/2, - h/2), in Eq. Following the 
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same steps, it is easy to arrive at an equivalent approxi- 
mation for the other missing pair, G{t + h/2, + h/2) + 
G{t — h/2,j\ + h/2). We summarize these below: 

G{t + h/2, - h/2) + Git - h/2, n - h/2) = 

-Vin-h/2) {f{t,r,-h) + f{t,r^)) 



G{t + h/2, n + h/2) + G{t - h/2, n + h/2) = 

- V{r, + h/2) {fit, + h) + fit, n)) 



2 V 2 



1 fh 



+ 



'-2^2 
Vir^ - h/2) 



Vir^ - h/2) 

S„!i dudv + Oih'^) 



Vir, + h/2) J J 



5,,, dudv + Oih'^). 

(31) 



Clott 



(30) 



Except for the presence of integrated source terms, these 
equations are identical to Lousto's equations (32) and 
(33), and Haas's equations (2.8) and (2.9). 



J 



Following Haas [Tj], we choose to avoid isolated occurrences of /(t, r*), which prove to be numerically unstable 
close to the event horizon. As pointed out in [ij], this is due to having first approximated G = —Vf, which makes 
it difficult to isolate / = —G/V where V ^ 0. This appears unnecessary if / were directly approximated instead of 
G in ((23)) . Nevertheless, like Haas, we avoid needing to isolate / by adding up equations Eq. pO]) and Eq. (ISTI) . and 
then Taylor-expanding the potential terms that are multiplied by /(f, r*). The result is Haas's equation (2.10) with 
extra source terms: 

^ G = G(t + h/2,r^ - h/2) + Git - h/2, - h/2) + Git + h/2, + h/2) + Git - h/2, n + h/2) = 



2Vin)fit,r,) 



1 



1 /h 



2 V 2 



Vir.,~h/2)fit,r.^-h) 



1 fh 



2\2 



Vir^ - h/2.) 



Vir,+h/2)fit, r, + h) 
1 



Vir.,+h/2) 



[Vir, - h/2) - 2Vir,) + Vir, + h/2)]ifit, - h) + fit, r, + h)) 



y(n - h/2) 



S^tt du dv - 



y(r, h/2) 



Cloft 



S,,, dudv + Oih'^) 



(32) 



Cright 



This last equation completes the pieces needed for the evolution algorithm. 

Using Eq. (flQ)) for the derivative term and Eq. (|20|) for the potential and source terms, the result of integrating the 
wave equation over the unit cell finally yields: 

/(. ,h,r.).- fit -h,r.), [^r^y;±^P] /(., r.,h), 7 ^^I^ ' (^^^ /(., . ~ .) 



[i + idmr-*)] 



[i + airvir^)] 



) (4Go + dudv 



Oih% 



(33) 



where Go is evaluated according to Eq. with 
G(t,r*) = — y(r*)/(t, r,); ^ G is the expression in 
Eq. ([5^ ; and the double Simpson rule Eq. pp)) is applied 
in evaluating the remaining integral term JJ^ S^ttdudv. 
With this equation, one can now determine the field / at 
time t + h given its values at earlier times t and t ~ h. 

This derivation makes liberal use of double Simpson 
and double trapezoidal formulas when approximating in- 
tegrals of the source and potential terms over the unit 
cell. The formulas come from their single-integral coun- 



terparts: 

fi^)dx = - ifixo) + fixo + h)] iO 

(34) 

/ fix)dx = J [/(xo) -f 4/(xo + h) + fixo + 2h)] 

Jxo 

-'^f^'HO. (35) 
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where /*^"^ denotes the nth-derivative of /, and ^ is some 
point within the hmits of integration. These require the 
boundedness, if not existence of the second and fourth 
derivatives of the integrand for the error estimate to be 
vahd. With the limited diffcrentiabihty of our source 
and potential terms (C° and C^, respectively), one might 
worry about the validity of our over-all convergence es- 
timate. However, our calculations reveal that 4th-order 
convergence is achieved despite this deficiency. 



D. Initial data and boundary conditions 

For our evolution we have no obvious method for choos- 
ing a priori the correct initial data, which consists of the 
value of "0^ on two consecutive constant-time slices. Con- 
sequently we just set the initial i/j^ fo zero everywhere 
on the initial two slices. Physically, this scenario corre- 
sponds to the impulsive appearance of the scalar point 
charge along with , which leads to spurious radia- 

tion contaminating our computational domain during the 
early stages of the evolution. Fortunately, this radiation 
propagates out of the regions of interest quickly; so to 
circumvent the need for proper initial data, we simply 
evolve the equation to long enough times such that ini- 
tial data effects do not become pertinent in any of our 
results. 

With the scalar charge moving in a circular orbit, it 
is expected that the field eventually becomes stationary 
in a frame corotating with the charge. A practical test 
then for the persistence of initial data effects is to simply 
check whether or not the field has already settled into a 
quiescent state when evaluated in this frame. 

Boundary conditions are treated similarly. Rather 
than handling them carefully, we instead made the com- 
putational domain large enough that errors incurred by 
unspecified boundary conditions did not affect our re- 
gions of interest. For this work, our choice of boundaries 
were at = -700M and = 800M. 



E. Self-force calculation 

At the end of evolution for each mode, we compute the 
self-force at the location of the particle. Since this loca- 
tion is not on any grid point, interpolation of /;„i(T, r) 
and its derivatives to r = i? was required using a selection 
of grid points surrounding it. 

Once this was done, computing the self-force was a 
simple matter of performing the following sums: 



1=0 



L I 



R 



1=0 m=l 

(37) 



1=0 

L I 

+ E E I^'^ {9rflm{T, R)Yi„, (|, 



1=1 m=l 



(38) 



Here, L is the point where we truncate the multipole ex- 
pansion. In all our work we have used L — 39. These 
sums arise primarily because our charge moves in a cir- 
cular orbit. 

Two methods were employed for interpolation. The 
first was a simple Lagrange interpolation of both / and 
drf to r = R. However, because of the finite differentia- 
bility of our regular field at r = i? we also interpolated 
using the form 

V'^(r) =Ao + Aix + A2x^ + A^x^ + 9{x)Bqx^, (39) 

where x = r — R, and 9{x) is the standard Heaviside 
function. This form closely respects the nature of the 
regular field at r = i? by allowing for a discontinuity in 
the third derivative. 

With this form, (Vr-F)|r=i? — Ai. However, this led to 
results not significantly different from the one achieved 
with ordinary Lagrange interpolation. 



V. CODE DIAGNOSTICS 
A. Convergence 

The convergence of a time-domain code is easily deter- 
mined by computing the convergence factor n as defined 
by Loustoflil: 



n{r^, t) =log 



fAh{r*-t) - f2hir*,t) 



f2h{r*,t) - fh{r*,t) 
-Hog|e(")(C)|/log(2), 



/log(2) 



(40) 



1=1 m=l 



where /a (r, , t) is the result of the evolution for a res- 
olution of A, and e'-"-'(C) represents an error function 
« 1. An nth-order evolution code is one for which 
V' = ^JN{h) + (e^"'')(^)/i" , where ipN{h) is the numer- 
ical solution at resolution h. 

In checking convergence, one evolves the wave equation 
at different resolutions, h, 2h, and 4ft,. For a fixed — R, 
one then extracts fh{R,t), f2h{R,t), and f4h{R,t) for all 
t. From these, one can compute n{R,t). 

The convergence factor was computed for a few repre- 
sentative points in the wavezone and in the region close 
to the point particle. Two of these are shown in Fig. [51 



9 



These are for r « lOM and r « lOOM. All show the de- 
sired 4th-order convergence eventually, following a tran- 
sient period in which the numerical evolution is contam- 
inated by the effects of poor initial data. 



convergence at r=10M 
convergence at r=100M 



100 



200 



300 



400 
t/M 



500 



600 



700 



FIG. 5: Convergence at the particle location (r = lOA/) and 
in the wavezone (r = lOOM). At the start of the evolution, 
inequivalent initial data lead to the lack of 4th-order conver- 
gence. But n gradually approaches 4 as initial-data effects 
propagate away from the computational domain. Note that 
the convergence test at the particle location already includes 
the interpolation step. 



B. High-Z fall-off 

In a i a i, 

it was demonstrated that the rate of 
convergence of the /-components of the self-force was dic- 
tated primarily by the lack of differentiability of the regu- 
lar piece from which the self-force is computed. By defini- 
tion, the difference between the retarded field and the sin- 
gular field yields a function that is C°° . In this ideal sit- 
uation, convergence in I of the self-force computed from 
this smooth regular field would be exponentially fast. In 
practice, however, one is always limited to constructing 
only an approximate singular field, therefore leaving non- 
differentiable pieces in the residual V'™* — V''^- The degree 
of non-differentiability of this remainder is what sets the 
rate of convergence of the self- force in I. 

The high-/ asymptotic structure of the singular piece 
tp^ is such that: 



lim(V^V )i = [l + -]Ar + Br 



(2Z- l)(2/-f 3) 



(2/-3)(2/- l)(2Z-h3)(2Z-H5) 



(41) 



where A^B,D,... are the regularization •parameters, 
which commonly appear in contemporary self-force stud- 
ies 

The number of regularization parameters that can be 
determined in this expansion corresponds directly to the 
accuracy of the singular field approximation. Conver- 
gence in / of the self-force V^V'^ is then fixed by the 
lowest-order undetermined piece of the approximate sin- 
gular field. Specifically, if the singular field is accu- 
rately determined only up to the i?-term of the expansion 
above, then the /-convergence of the self-force would be 
^ 1//^, corresponding to the D-term fall-off. 

The approximation to the singular field here is '0'^ = 
q/p, and the attendant THZ-Schwarzschild coordinate 
transformation, has been shown in Q to include at least 
the _D-term. The expectation then would be for the l- 
components of our remainder, (VrV'^);; to off ^is the 
piece: 



/■^3/2 



(2/-3)(2/- l)(2/-h3)(2/ + 5)' 



(42) 



Fig.[6]shows our results confirming this expectation. Our 
results are plotted with L>, and S^^) fall-off curves 
found in Eq. (|4T|) that arc made to match at / = 15. 
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FIG. 6: (VrV'^)i versus /. Our results show /-convergence 
closest to the _E'^' fall-off. The blue lines correspond to the 
expected fall-off in the r-component of the self- force when one 
regularizes using a singular field approximation without the 
D-term, iJ'-^-'-term, and .E'-^-'-term, respectively. Our result is 
matched to these curves at /=15. 

The t-component of the self-force, on the other hand, 
does not require regularization for the case of a charge in 
a circular orbit of Schwarzschild. An exponential fall-off 
is then expected. This is shown in Fig. [T] 
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FIG. 7: {Vt'ip^)i versus I. The expected exponential fall-off 
is observed until the point where numerical noise begins to 
dominate. 



C. Dependence on the window function 

The use of a window function 1^ is a peculiar feature 
of our approach. Its function is mainly to kill off ip^ in 
the regions where it is no longer relevant and thereby to 
have the computed regular field "0^ transform into the 
retarded field in those regions. As it is a mere artifact of 
our implementation, it is crucial that the self-force and 
waveform be independent of the specific choice of window 
function. 

One has considerable freedom in choosing W , the only 
requirements being that W goes to 1 and that its gradient 
vanishes fast enough in the limit that one approaches the 
point particle. With our specific choice of W becoming 
numerically significant only in an annular region |r— i?| < 
tJ, we have inspected the changes in the self-force and 
fluxes as one varies the width a. 

Fig. [5] shows the effect of doubling the annular sup- 
port of the window function. The (/ = 2, m = 2) wave 
equation was evolved for the same length of time, but 
with effective sources having different window functions. 
A comparison is then made of the resulting fields over 
most of the computational domain. It is seen that the 
fields differ significantly only in regions where the win- 
dow functions differ. Nevertheless, the regular field ip^ 
remains the same (up to fractional changes of ~ 10~® in 
the most physically- relevant regions: the vicinity of the 
charge, r = lOM (where the self- force is computed), and 
the wavezone, r ^ lOM (where the waveform is to be 
extracted) . 

As desired then, the window function appears to have 
no effect on any of the numerical results attained. 




-20 -10 10 20 30 40 50 60 70 80 90 100 110 12 
r./M 



FIG. 8: Fractional changes in ip22 as a result of using different 
window functions. Note that these changes are significant 
only where the window functions differ; they are insignificant 
in the important regions in the vicinity of the charge, r = 
lOM, and in the wave zone, r 3> lOM. 



VI. RESULTS 

A. Recovering the retarded field 

From our numerical calculations we are able to accu- 
rately recover the retarded field. In the wavezone, where 
the singular field is negligible, this retarded field equals 
our regular field ij:^. Since, energy fluxes depend directly 
on the retarded field in this region, the accuracy with 
which we recover the retarded field in the wavezone gives 
us a measure of how well we can compute fluxes using our 
method. We determine this accuracy by comparing our 
result for t/;^ in the wavezone with that obtained for the 
retarded field using a separate frequency-domain calcula- 
tion. An example of such a comparison is shown in Figs. [5] 
andfTOl We observe relative errors that are at worst 10~^. 
Shown in Fig. [H] are the (/ = 2, to = 2) component of the 
retarded field computed in the frequency-domain and our 
corresponding time-domain result, ^22 + {Wip^)22, for the 
case of a charge at r = lOM . Also shown are the singular 
field {W 11)^)22 and the regular field ^/i^. 



B. Self- force 

We obtain the t and r components of the self-force 
for an orbit at radii R = lOM and 12M. These are 
summarized in Table HI 

Fig. [TT] shows the convergence of our calculation of the 
the {I = 2,m = 2) time component of the self-force we 
show the convergence of our time-domain calculation to 
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FIG. 9: Comparison of time-domain and frequency-domain 
results for /22 (''»)• The regular field is the result of our code 
(represented by the blue dashed line). Adding this to the 
(?=2,m=2)-component of our analytical singular field, W-t/)^ , 
results in the FD-computed retarded field to good agreement. 
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FIG. 10: Relative error between time-domain and frequency- 
domain results for /22(r*). Excellent agreement is achieved; 
errors are at worst ~ 10~®. 



the frequency-domain result. We have excellent conver- 
gence after a time of 200M, which is approximately one 
orbital period. 





R 


Time-domain 


Frequency-domain 


error 




lOM 3.750211 X 10"" 


3.750227 X 10" 


-b 


0.000431% 




lOM 


1.380612 X 10~" 


1.378448 X 10-^ 


-5 


0.157% 




12M 


1.747278 X 10"'' 


1.747254 X 10" 


'5 


0.00139% 




12M 


5.715982 X lO"*^ 


5.710205 X 10" 


-6 


0.101% 



TABLE L Summary of self-force results for R = lOM and 
R = 12Af . The error is determined by a comparison with an 
accurate frequency-domain calculation [^. 
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FIG. 11: Relative error in the time-domain calculation of /22, 
as compared with the frequency-domain calculation, versus 
time at r close to the charge at lOM. 



VII. DISCUSSION 

In the specific context of a point charge orbiting a back 
hole, we have introduced a very general approach suitable 
for the time-domain generation of waveforms and also the 
calculation of the backreacting self-force. 

Our initial tests are admittedly on the very restric- 
tive case of circular orbits of the Schwarzschild geometry, 
where we have taken advantage of the spherical symme- 
try to decompose the source and field into spherical har- 
monics. This has allowed us to compare our results to 
available frequency-domain results of very high precision 
[l,[l2|- In some manner, because of our use of a spherical- 
harmonic decomposition, our analysis might be likened to 
using spectral methods. But, the method of field regular- 
ization inherently does not require a mode-decomposition 
and could be implemented with a full (3-1-1) numerical 
code. For our test case we achieve an extremely accurate 
calculation for the time component of the self-force dtip^, 
which is equivalent to the rate of energy lost by radia- 
tion. Notably, in our (1+1) implementation, the initial 
data settled down to provide this accurate component 
of the self-force within only one orbit of the particle as 
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shown in Fig.[TTJ This might be contrasted with a calcu- 
lation of dE/dt made from a flux integral evaluated in the 
wave-zone, which, with similarly unspecified initial data, 
requires evolution over a substantial number of orbits. 

For a circular orbit, the radial component of the self- 
force drip^ is conservative and generally more difficult to 
calculate. We were able to match more accurate anal- 
yses i, El, [H [11 to about 0.1%. With the spher- 
ical harmonic decomposition, our analysis went up to 
L = 39. This relatively high number is due primarily 
to the slow polynomial convergence resulting from the 
mode-decomposition of the self-force. We expect this to 
be endemic in all self-force calculations that rely on some 
kind of spectral decomposition, as it is the penalty in- 
curred when one represents objects of limited differentia- 
bility in terms of smooth functions. 

A technique similar to that described in Ref. could 
possibly mitigate this weakness. For our specific imple- 
mentation, we could choose to calculate and sum modes 
only up to, say / = 15, and then take advantage of the 
known asymptotic fall-off in I shown in Eq. (j4ip . Using 
the computed modes, we determine the coefficients in the 
expected fall-off for the self-force, and then using these, 
analytically complete the sum to Z = oo. This results in 
a slightly more accurate result for drip^- A similar pro- 
cedure of "fitting" to a known asymptotic fall-off might 
prove useful if one chooses to implement field regulariza- 
tion using spectral methods. 

We expect field regularization to be best implemented 
on a (3-1-1) finite-difference code, with mesh refinement 
in the vicinity of the charge to better resolve the limited 
differentiability of our analytically constructed source 
function. Such a process will ameliorate the problem 
of slow polynomial convergence ailing typical mode-sum 
prescriptions. 

For the EMRI problem today, there is great interest in 
calculating the rate of energy being radiated for a point 
mass orbiting a rotating black hole and in using the result 
to modify the the orbit of the mass with some version of 
an adiabatic approximation. For a general orbit, the en- 
ergy flux is not easy to determine. Current methods use 
the axial symmetry of the Kerr geometry to separate out 
one dimension, and then deal with a (2-1- 1)D problem for 
the radiation from a point mass. The representation of a 
point mass on a grid is typically problematical. Replac- 
ing a (5-function source by a narrow Gaussian [2^, [28| is 
reasonable but docs not accurately reproduce frequency- 
domain results. The recent distribution of a ^-function 
over a modest number of grid points by Sundararajan 
et al t23j appears more robust. The strategy laid out 
here provides a natural remedy to this issue. Instead of 
dealing with a wave equation with a ^-function source, 
we solve an equivalent problem with a regular and dis- 
tributed source. The results displayed in Figures ^ and 
[TU] clearly point to the effectiveness and accuracy of our 
method. 

Beyond the numerical modeling of (5-function sources 
though, the method of field regularization provides di- 



rect access to the self-force, which is essential in a fully- 
consistent treatment of particle motion and wave gen- 
eration. Current methods under development are based 
upon energy and angular momentum flux calculations 
that will certainly miss conservative self-force effects. 
These methods rely upon flux integrals evaluated in the 
wave-zone and some orbit averaging or post-processing to 
effect the change in orbital energy or angular momentum, 
which are difficult to implement carefully [13, [IJ and to 
justify rigorously. The more direct approach of locally 
calculating the self-force to update the particle orbit has 
been largely avoided because of the prohibitive compu- 
tational expense associated with mode-sum calculations 
of the self-force. In a (3-1-1) finite-differencing implemen- 
tation of field regularization, calculating the self-force is 
no more expensive than performing a numerical deriva- 
tive and possibly an interpolation. As such, it represents 
a step forward towards the goal of efflcicntly producing 
consistent numerical models of particle motion and radi- 
ation in curved spacetime. 

The recent proposal of Barack, Golbourn and Sago 
[isl . [H, [l^l is closest in spirit to our method of field regu- 
larization. They model a point charge with a distributed 
effective source derived instead from their "puncture 
function", which is quite similar to our tp^ ■ They base 
their construction of the puncture function on the 'di- 
rcct'-|-'tair decomposition, rather than on the Green 
function decomposition in [l9| that naturally provides 
our regularizing singular field . Their current punc- 
ture function, however, appears to prevent them from 
calculating a self-force. Moreover, in anticipation of 
a Kerr background application, they envisage using a 
(2-1-1) code, necessitating a mode-sum over a mode in- 
dex TO, which will again feature the characteristic poly- 
nomial convergence of this approach to self-force calcula- 
tion. This is demonstrated in Fig. [T^l Using our results, 
we perform partial sums over I, i.e. 

L 

(V,V^«)„ EE {^ri^%m, (43) 

i>|m| 

to get the resulting fall-off in m. We observe a fall-off 
close to l/m^ in the modes. Consequently, if we were to 
follow Barack et al's m-mode prescription, the self-force 
would converge as 1/m'^. 

In principle, our method of field regularization ap- 
pears to resolve two important issues in the context 
of EMRI simulations: (a) numerically representing i5- 
fimction sources, and (b) calculating the self-force. At 
this time our test of a scalar charge in a circular orbit 
of the Schwarzschild geometry is a carefully controlled 
numerical experiment and provides us with detailed in- 
formation about the relationship between the approxi- 
mation for tp^ and the rate of convergence of the self- 
force. However, our test is also extremely elementary 
when compared to the actual case of a point mass emit- 
ting gravitational waves from a generic orbit of the Kerr 
geometry, which is most relevant for EMRIs. Future work 
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FIG. 12: {Vrip^)m versus m. Our results show an m-fall-off 
closest to 1/rn^. 



will focus on exploring the robustness of our technique 
against these more interesting cases. 
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APPENDIX A: APPROXIMATE SINGULAR 
FIELD IN THZ-COORDINATES 

A special coordinate system developed by Thorne and 
Hartle |2fl(| and by Zhang [2l| is particularly useful for 
self-force analyses. These THZ coordinates {t,x,y,z) 
are defined in a neighborhood of a geodesic of a vac- 
uum spacetime and are locally incrtial, harmonic and 
Minkowskii-like, and centered on the geodesic with t mea- 
suring the proper time along the geodesic. In these spe- 
cial coordinates our expression for the approximate sin- 
gular field ^ps — q/ yj \ + is quite simple. But, 
for the case of a point charge in a circular orbit about a 
Schwarzschild black hole this simplicity of t/i"^ belies the 
hidden complexity of the coordinate transformation be- 
tween the Schwarzschild coordinates (is,r, 0, (/)) and the 
THZ coordinates. 

The full coordinate transformation, which we use for 
our analysis in the main body of this paper may be found 
in Eqs. (B1)-(B9) of Ref. Below we give only an 
abbreviated form of these equations to give a sense of 
how the coordinate transformation is implemented. The 
formulae below give x, y, z and t as smooth functions of 
the Schwarzschild r, 6*, </> and t^. We can then define a 
function p — + y'^ + which has the property that 



VVail/p) = -^n5{x) + 0[l/p). 



(Al) 



If we had used q/ p as the approximation -08 for the sin- 
gular field then the effective source for the regular field 
-0^ in the vicinity of the point charge would be singular. 



off ■ 



-V''Vb{q/p)-A7rqSix)^0{l/p), 



(A2) 



rather than 0{p), which is the case for the effective source 
which we actually use as described in Eq. PU)) . 

The coordinates which lead to p are now given for 
a circular geodesic of the Schwarzschild geometry at 
Schwarzschild radius R: We first define two useful func- 
tions 



[r sin 6 cos((/) — fits) — -R] 



M 



(1 - 2M/i?)i/2 
(r - i?)2 



[2(1 



2M/R) 
+ 0{p') 



R^{l-2M/Ry/^ 
R^ sir? e sir? [4) - fits) + R^ cos^ 6 

(A3) 



and 



y 



r smf sm 



R~2M 
R~3M 



1/2 

+ 0(p3).(A4) 



The 0{p^) terms indicate that these (and the formulae 
below) could be modified by the addition of arbitrary 
0{p^) terms without necessarily changing the usefulness 
of these coordinates. 

In terms of these two functions, the THZ coordinates 
(i, X, y, z) are 



x cos 



y sm 



(A5) 
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and 



y — X s\n{VlHs) + y cos(f2^ts) 



where Vt'^ = Vl^Jl - 3M/R, along with 
z = rcos{9) +0{p^) 



(A6) 



(A7) 



and 



t = 



rilR sin 9 sin((/) — fiig) 
R-3M 



+ 0(p3) (A8) 



The set of functions i, y, z) forms a non-inertial co- 
ordinate system that co-rotates with the particle in the 
sense that the x axis always lines up the center of the 
black hole and the center of the particle, the y axis is 
always tangent to the spatially circular orbit, and the z 
axis is always orthogonal to the orbital plane. 

The THZ coordinates (t, x, y, z) are locally inertial and 
non-rotating in the vicinity of the charge, but these same 
coordinates appear to be rotating when viewed far from 
the charge as a consequence of Thomas precession as re- 
vealed in the dependence in Eqs. (jASp and (|A6[) 
above. 

The coordinates {t,x,y,z) given above are said to be 
second order THZ coordinates and differ from the actual 
fourth order ones used in the main body of this paper by 
the replacement of the 0{p^) terms appearing above by 



specific terms which scale as and @ and leave the 
undetermined parts of the THZ coordinates being 0{p^). 



APPENDIX B: CONVERGENCE FACTOR 

An nth-order evolution code is one for which — 
ij}N{h) -t- (e(")(^))ft," , where ipNih) is the numerical so- 
lution at resolution /i, and e*-"-* (^) is some unknown error 
function or order w 1. Consider three resolutions /i, 2h, 
Ah. This then leads to 



Thus, 



and so 



\^l;N{2h)-^M{h)\ ' 



(C)|2", 



(Bl) 
(B2) 
(B3) 



(B4) 



log 



V'7v(4/i)-Vjv(2/i) 



^N{2h)-^N{h) 
-log|e(")(0|/log(2) 



/log(2) 



(B5) 
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